#-----------------------------------------------------------+
#                                                           |
# do_all_fisher.R - Do all Fisher tests for occurrence data |
#                                                           |
#-----------------------------------------------------------+
#                                                           |
#  AUTHOR: James C. Estill                                  |
# CONTACT: JamesEstill_at_gmail.com                         |
# STARTED: 11/12/2008                                       |
# UPDATED: 11/14/2008                                       |
#                                                           |
# DESCRIPTION:                                              |
#   Do all fisher tests for all paired sets of column 2     |
#   data given input in text file as:                       |
#                                                           |
#     column_1 <tab> column_2                               |
#                                                           |
#   where column_1 is data such as BAC_ID and column_2 is   |
#   data such as TE_occurrence_id.                          |
#   This will test for significant positive or significant  |
#   negative associations of the col_2 data on the col_1    |
#   objects.                                                |
#                                                           |
#-----------------------------------------------------------+

#-----------------------------+
# CLEAR WORKSPACE             |
#-----------------------------+
rm(list = ls(all = TRUE));

#-----------------------------+
# OPEN INPUT FILE             |
#-----------------------------+
#infile <- "/Users/jestill/projects/gina_count/test_in.dat";
infile <- "/Users/jestill/projects/gina_count/bac_nolowcopy.dat";
#infile <- "/Users/jestill/projects/gina_count/BAC_Id_and_fam.dat";
#infile <- "/Users/jestill/projects/gina_count/all_bac_data.dat";

in_data <- read.table(infile, header=F);

#-----------------------------+
# OUTPUT FILE PATHS & OPTIONS |
#-----------------------------+

# SUMMARY OUTPUT FILE
out_file <- "/Users/jestill/projects/gina_count/no_low_associations.txt";
# Write output file headers
cat (c ("Val_1\tVal_2\tDirection\tP_value\tEstimate\tNullVal\tBF_p_val\n"  ),
     file=out_file, append=FALSE);
# Print fisher tables to output file (SET TO FALSE OR TRUE)
print_fisher = TRUE;

# P VALUE MATRIX (UNCORRECTED)
p_mat_out_file <- "/Users/jestill/projects/gina_count/no_low_bac_p_vals.txt";

# ESTIMATED ODDS RATIO
est_mat_out_file <- "/Users/jestill/projects/gina_count/no_low_bac_est_vals.txt";

#-----------------------------+
# GENERATE CONTINGENCY TABLE  |
#-----------------------------+
# Summarize all data as contingecy table
cont_table <- table(in_data);

# Can also get the table dimensions
tab_dim<-dim(cont_table);
num_rows <- tab_dim[1];
num_cols <- tab_dim[2];

# Will need to get the names of the rows and the cols
# This will be used to iterate across the list
tab_row_names <- rownames(cont_table);
tab_col_names <- colnames(cont_table);

# THE NUMBER OF TESTS IS EQUAL TO THE TOP TRIANGLE OF THE COMPLETE MATRIX
# SELF COMPARISONS NOT INCLUDED
# This is used for bonferroni corrections below
num_tests <- ( ((num_cols^2) - num_cols)/2 );

#-----------------------------+
# INITIALIZE MATRICES         |
#-----------------------------+

# P VALUE MATRIX
p_val_mat <- matrix(data=0,nrow=num_cols, ncol=num_cols,
                    dimnames= list(tab_col_names,tab_col_names) );

# ESTIMATED ODDS RATIO MATRIX
est_val_mat <- matrix(data=0, nrow=num_cols, ncol=num_cols,
                      dimnames= list(tab_col_names,tab_col_names) );

# NULL VALUE MATRIX
# This is used when generating the info for multiple comparisions at the end
null_val_mat <- matrix(data=0, nrow=num_cols, ncol=num_cols,
                      dimnames= list(tab_col_names,tab_col_names) );

# FISHER TABLE ARRAY
# 3-D array to hold all of the fisher table counts
if (print_fisher == TRUE) {
  fisher_table_ary <- array(data=NA, dim=c(num_cols,num_cols,4) );
}

# P VALUE VECTOR
# This will be used for sequential bonferroni
p_val_vect <- vector( mode="numeric", length=num_tests);
p_count <- 0; # To increment for proper place in p_vector

#-----------------------------+
# DO ALL FISHER TESTS         |
#-----------------------------+

# Make fisher result matrix

for (i in 1:num_cols) {
  for (j in 1:num_cols) {

    fish_table <- matrix(data=0,nrow=2,ncol=2);


    # FISHER TABLE AS:
    #               j_present | j_absent
    #            +------------+-------------+
    #   i_present|     1      |      2      |
    #            |   [1,1]    |    [1,2]    |
    #            +------------+-------------+
    #   i_absent |     3      |      4      |
    #            |   [2,1]    |    [2,2]    |
    #            +------------+-------------+
    # fish_table_ary 3rd dim position indicated above in cells
    
    for (n in 1:num_rows) {

      #-----------------------------+
      # FILL FISHER TABLE           |
      #-----------------------------+
      # Doing evaluation as  > 0 and not = 1 takes into account cases where
      # there is more the one copy of col 2 data on col 1 data
      # For example, there may be more then one copy of an element
      # on a BAC and you just want to treat this a presence absence.

      if (cont_table[n,i] > 0) {
        if (cont_table [n,j] > 0) {
          # increment 1,1 (+|+)
          fish_table[1,1] = fish_table[1,1] + 1;
        }
        else {
          # increment 2,1 (+|-)
          fish_table[2,1] = fish_table[2,1] + 1;
        }
      }
      else {
        
        if (cont_table [n,j] > 0) {
          # increment 1,2 (-|+)
          fish_table[1,2] = fish_table[1,2] + 1;
        }
        else {
          # increment 2,2, (-|-)
          fish_table[2,2] = fish_table[2,2] + 1;
        }
        
      }
      
    } # End for for each row

    #-----------------------------+
    # DO FISHER'S EXACT TEST      |
    #-----------------------------+
    fish_result <- fisher.test(fish_table, conf.int=TRUE);

    # Possible results to fetch from Fisher's Test
    #fish_conf <- fish_result$conf.int;     # confidence interval for odds ratio
    #fish_p_val <- fish_result$p.value;     # P value
    #fish_estimate <- fish_result$estimate; # Estimate of odds ratio
    #null_odds <- fish_result$null.value;   #

    # LOAD RESULTS TO MATRIX
    est_val_mat[i,j] <- fish_result$estimate;
    p_val_mat[i,j] <- fish_result$p.value;
    null_val_mat[i,j] <- fish_result$null.value;

    # LOAD FISHER TABLE COUNTS TO ARRAY
    if (print_fisher == TRUE) {
      fisher_table_ary [i,j,1] <-fish_table[1,1];
      fisher_table_ary [i,j,2] <-fish_table[1,2];
      fisher_table_ary [i,j,3] <-fish_table[2,1];
      fisher_table_ary [i,j,4] <-fish_table[2,2];
    }

    #-----------------------------+
    # REPORT SIGNIFICANT          |
    # ASSOCIATIONS                |
    #-----------------------------+

    # Only report for i < j
    if (i < j) {

      # LOAD P VALUE TO VECTOR
      p_count <- p_count + 1;
      p_val_vect[p_count] <- fish_result$p.value;
    
      if (fish_result$p.value < 0.05) {

        # Bonferroni corrected for multiple comparisons
        bf_corrected_p <- ( num_tests * fish_result$p.value );
          
        # This can be modified with a different null then the default
        # and the code will still work.
        if ( est_val_mat[i,j] > null_val_mat[i,j] ) {
          cat (c (tab_col_names[i],"\t",tab_col_names[j],"Positive\t",
                  fish_result$p.value,"\t",fish_result$estimate,
                  "\t",fish_result$null.value,
                  "\t",bf_corrected_p,"\n"  ),
               file=out_file, append=TRUE);
          
          # The contignecy table
          if (print_fisher == TRUE) {
            cat (c (fish_table[1,1],"\t", fish_table[1,2], "\n",
                    fish_table[2,1], "\t",fish_table[2,2], "\n" ),
                 file=out_file, append=TRUE);
          }
        }
        else {

          # Print Fisher Test Results to outfile
          cat (c (tab_col_names[i],"\t",tab_col_names[j],"Negative\t",
                  fish_result$p.value,"\t",fish_result$estimate,
                  "\t",fish_result$null.value,
                  "\t", bf_corrected_p,"\n" ),
               file=out_file, append=TRUE);
          
          # The contignecy table
          if (print_fisher == TRUE) {
            cat (c (fish_table[1,1],"\t", fish_table[1,2], "\n",
                    fish_table[2,1], "\t",fish_table[2,2], "\n" ),
                 file=out_file, append=TRUE);
          }
          
        } 
        
      } # End if p value < 0.05
    } # End of if i<j
    
  } # End of iterate across j
  
  # Stop, just do first 
  #stop ("Just testing");
} # End of iterate across i


#-----------------------------+
# WRITE MATRIX OUTPUT FILES   |
#-----------------------------+

# Matrix of the estimates of odds ratios
write.table (est_val_mat, file=est_mat_out_file);

# Matrix of P values
write.table (p_val_mat, file=p_mat_out_file);

#-----------------------------+
# ALTERNATIVES TO BONFERRONI  |
#-----------------------------+
# If the p values above are loaded into a vector
# it is possible to use the p.adjust to adjust
# p values
# see help("p.adjust") for methods to adjust the p value
p_adjusted <- p.adjust (p_val_vect, method="holm");

# PRINT HEADER TO SHOW THAT WE ARE SHOWING HOLM RESULTS
cat ( c("\n\n +------------------------------------------------+\n",
        "| MULTIPLE COMPARISON CORRECTED                  |\n",
        "+------------------------------------------------+\n\n"
        ),
     file=out_file, append=TRUE);

cat (c ("Val_1\tVal_2\tDirection\tP_value\tEstimate\tNullVal\tBF_p_val\tHolm_p_val\n"  ),
     file=out_file, append=TRUE);


#-----------------------------+
# REPORT SIGNIFICANT          |
# ASSOCIATIONS                |
#-----------------------------+
# The following will only report p-values for paired sets that
# are significant when adjusted for multiple comparisons.
p_count <- 0;
for (i in 1:num_cols) {
  for (j in 1:num_cols) {

    if (i < j) {
      
      p_count <- p_count+1;
      
      if (p_val_vect[p_count] < 0.05) {
        
      # Bonferroni corrected for multiple comparisons
        bf_corrected_p <- ( num_tests * p_val_mat[i,j] );
      
        if ( est_val_mat[i,j] > null_val_mat[i,j] ) {
          cat (c (tab_col_names[i],"\t",tab_col_names[j],"Positive\t",
                  p_val_mat[i,j],"\t",est_val_mat[i,j],
                  "\t",null_val_mat[i,j],
                  "\t",bf_corrected_p,
                  "\t",p_val_vect[p_count],"\n" ),
               file=out_file, append=TRUE);
          
          # The contignecy table
          if (print_fisher == TRUE) {
            cat (c ( fisher_table_ary [i,j,1], "\t", fisher_table_ary [i,j,2], "\n",
                     fisher_table_ary [i,j,3], "\t", fisher_table_ary [i,j,4], "\n" ),
                 file=out_file, append=TRUE);
          } # End if print_fisher is TRUE
        }
        else {

          # Print Fisher Test Results to outfile
          cat (c (tab_col_names[i],"\t",tab_col_names[j],"Negative\t",
                  p_val_mat[i,j],"\t",est_val_mat[i,j],
                  "\t",null_val_mat[i,j],
                  "\t", bf_corrected_p,"\t",
                  "\t", p_val_vect[p_count],"\n"),
               file=out_file, append=TRUE);

          if (print_fisher == TRUE) {
            cat (c ( fisher_table_ary [i,j,1], "\t", fisher_table_ary [i,j,2], "\n",
                    fisher_table_ary [i,j,3], "\t", fisher_table_ary [i,j,4], "\n" ),
                 file=out_file, append=TRUE);
          } # End if print_fisher is TRUE
          
        } 
        
      } # End if p value < 0.05
    } # End of if i<j
  } # End of iterate i
} # End of iterate j


#-----------------------------------------------------------+
# HISTORY                                                   |
#-----------------------------------------------------------+
#
# 11/12/2008
# - Base program written to input two col matrix,
#   transform to a contingency table and iterate
#   across columns to genereate Fisher tests.
#   Significant results of fisher tests are written
#   to an outfile.
# 11/14/2008
# - Cleaned up code
# - Added alternatives to straight bonferroni for the
#   multiple comparisons correction. This required
#   loading the p-values into a vector and using the
#   R p.adjust function.
# - Added print fisher tables to the multiple comparisions
#   corrected results. This required creating a 3D array to
#   hold all of the Fisher table counts.
